Lecture 5

Date: Sept 16, 2024

The goal of this lecture is to focus on inexact number types in Sage.

Inexact fields

These are described in the documentation. See for more details.

Number representations

Let $\beta \geq 2$ be an integer. Then we can represent each non-negative integer by a finite sequence of numbers in the set $\{0,1,2,\ldots,\beta-1\}$. Namely, if $(a_n,a_{n-1}, \ldots, a_0)$ is a finite sequence with each $a_i \in \{0,1,2,\ldots,\beta-1\}$, it represents the number $$\sum_{k=0}^n a_k \beta^k=a_n \beta^n + \ldots + a_2 \beta^2 + a_1 \beta + a_0.$$ Numbers are ordered like in decimal in order from most signficant to least significant. When writing numbers in every day life, we assume that the most significant number in the sequence (i.e., $a_n$) is non-zero. Often computers store numbers using sequences of a fixed length, and in this case the most significant number stored can be a zero.

There some important special cases:

• Base 10 representations are said to be decimal representations.
• Base 2 representations are said to be binary representations.
• Base 16 representations are said to be hexadecimal representations. Often hexadecimal is abbreviated as hex.

Examples The number $14$ has binary representation $1110$ since $14=8+4+2$.

In base 16, we represent the numbers $\{0,\ldots, 15\}$ by switching to letters once we run out of digits. So, we use $$\{0,1,2,3,4,5,6,7,8,9,a,b,c,d,e,f\}$$ with $a=10$, $b=11$, $c=12$, $d=13$, $e=14$ and $f=15$. So for example, $415$ would be represented as 19f since $$415 = 1 \cdot (16)^2 + 9 \cdot 16 + 15.$$

Hex and bin

bin(14)

'0b1110'

hex(415)

'0x19f'

Floating point numbers

Let $\beta$ be a base, i.e., an integer greater than or equal to two. Then, just like in the decimal system every real number can be expressed as an infinite sum of the form $$x = \pm \left(\sum_{k=-\infty}^0 a_k \beta^k\right) \times \beta^n,$$ where each $a_k \in \{0,1,\ldots, \beta-1\}$. We also assume that $a_0 \neq 0$ except in the case that $x=0$.

This representation of $x$ consists of a sign ($\pm 1$), the mantissa $\sum_{k=-\infty}^0 a_k \beta^k$ and the exponent $n$.

A simplified and not entirely correct description of a how a computer stores double precision floatting point numbers follows. Base $2$ is used, so data is stored in bits. The sign takes up one bit. Some number of bit say $K$ are used to store the mantissa, which then has the form $\sum_{k=-K}^0 a_k \beta^k$. Note that because $a_0=1$ (unless $x=0$), this bit does not need to be stored. The precision of a number is the number of bits known, which is $K+1$.

Typically some number of bits say $E$ are used to store the exponent which means that roughly the exponent $n$ satisfies $$-2^{E-1}+2 \leq n \leq 2^{E-1}-1.$$

Note that there are only $2^E-2$ values available for $n$, which leaves two values available for handling special values, including zero, infinity, and "not a number".

Generally, a calculation yielding a number whose exponent exceeds the allowable range might produce the value $0$, $+\infty$, $-\infty$, or even a value like NaN (Not a Number). In Sage and other mathematical systems, we have access to representations of integers of arbitrary size (limited by memory available), so for some types of numbers this is not something we have to worry about.

Real double field

In most programming languages, floating point arithmetic is done with doubles. This is also what most processors are optimized to do, so it is very fast. The real double field is RDF.

A double consists of $64$ bits, $52$ bits of which are devoting to storing the mantissa, one bit for the sign, and the remaining $11$ bits are devoted to the exponent. The precision is $53$.

2^10

1024

RDF

Real Double Field

5.45.parent()

Real Field with 53 bits of precision

The above field is not RDF. To construct 5.45 as a double we need to do this:

RDF(5.45)

5.45

5.45.parent() == RDF

False

RDF.precision()

53

You can convert a number into RDF using RDF(number). For examples:

x = RDF(pi)
x

3.141592653589793

We can find out exactly what number is stored using the sign_mantissa_exponent method. This returns a triple (s, m, e) where s is a sign ($\pm 1$), the mantissa m in this case is an integer (whose binary representation has a number of bits equal to the precision), and the exponent e is an integer. The number represented is $$s \cdot m \cdot 2^e.$$

s,m,e = x.sign_mantissa_exponent()
s,m,e

(1, 7074237752028440, -51)

The quantity $s \cdot m \cdot 2^e$ is a rational number:

xx = s*m*2^e
xx

884279719003555/281474976710656

If we want to observe that it is an approximation of $\pi$, we can compute the numerical approximation:

xx.n()

3.14159265358979

Limits on number storage coming from the mantissa.

Because of the limits on the mantissa, the number $1 + 2^{-52}$ is the smalest number that is greater than one that can be stored:

x = RDF(1) + RDF(2^-52)
x

1.0000000000000002

You can see this is the smallest because the mantissa has the form $1000\ldots0001$:

bin( x.sign_mantissa_exponent()[1] )

'0b10000000000000000000000000000000000000000000000000001'

In particular, $1 + 2^{-53}$ can not be stored:

x = RDF(1) + RDF(2^-53)
x

1.0

x == 1

True

Limits on number storage coming from the exponent.

The exponent must satisfy $e \leq 2^{10} - 1$

x = RDF(2^(2^10 - 1))
x

8.98846567431158e+307

If you attempt to compute $2^{2^{10}}$ you get $+\infty$:

RDF(2^(2^10))

+infinity

Number can be stored with exponents down to about $-2^{10}<e$.

x = RDF(2^(-2^10-50))
x

5e-324

x.sign_mantissa_exponent()[1]

4503599627370496

Arbitrary Precision reals

Arbitrary precision reals are implemented with RealField. By default numbers have the same precision as doubles:

RealField()

Real Field with 53 bits of precision

This field is also represented as RR in Sage:

RR

Real Field with 53 bits of precision

5.45 in RR

True

5.45.parent() == RR

True

We can increase the prevision to $100$ bits, which is $100/\log_2 10$ digits.

F = RealField(100)
F

Real Field with 100 bits of precision

x = F(-pi)

x

-3.1415926535897932384626433833

Again you can use the sign_mantissa_exponent method:

s,m,e = x.sign_mantissa_exponent()
s,m,e

(-1, 995610453248924340922087778488, -98)

xx = s*m*2^e
xx

-124451306656115542615260972311/39614081257132168796771975168

xx.n()

-3.14159265358979

Here you can see the bits of the mantissa:

bin(m)

'0b1100100100001111110110101010001000100001011010001100001000110100110001001100011001100010100010111000'

And check that there are 100 of them. (The first "0b" does not count.)

len(bin(m))

102

Numbers with larger exponents can be stored in RR as compared to doubles (RDF):

F(2^(2^11))

3.2317006071311007300714876689e616

But there is still a limit:

show(F(2)^F(2^61))

\(\displaystyle 3.4280180247809630182481095324 \times 10^{694127911065419641}\)

show(F(2)^F(2^62))

\(\displaystyle +infinity\)

Here we compute $\frac{1000}{log_2 10}$ digits of $\pi$

RealField(1000)(pi)

3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127

Rounding modes

SageMath uses MPFR a floating point library. Most often, an exact solution of a calculation is not produced, in that case the software can only choose a number approximating the actual value. MPFR has multiple rounding modes that determine how this is done.

I want to point out the value in two in particular. You can choose to always round down or always round up. For examples:

RealField?

Docstring:     
   RealField(prec, sci_not, rnd):

   INPUT:

   * "prec" -- (integer) precision; default = 53 prec is the number of
     bits used to represent the mantissa of a floating-point number.
     The precision can be any integer between "mpfr_prec_min()" and
     "mpfr_prec_max()". In the current implementation,
     "mpfr_prec_min()" is equal to 2.

   * "sci_not" -- (default: "False") if "True", always display using
     scientific notation; if "False", display using scientific
     notation only for very large or very small numbers

   * "rnd" -- (string) the rounding mode:

     * "'RNDN'" -- (default) round to nearest (ties go to the even
       number): Knuth says this is the best choice to prevent
       "floating point drift"

     * "'RNDD'" -- round towards minus infinity

     * "'RNDZ'" -- round towards zero

     * "'RNDU'" -- round towards plus infinity

     * "'RNDA'" -- round away from zero

     * "'RNDF'" -- faithful rounding (currently experimental; not
       guaranteed correct for every operation)

     * for specialized applications, the rounding mode can also be
       given as an integer value of type "mpfr_rnd_t". However, the
       exact values are unspecified.

   EXAMPLES:

      sage: RealField(10)
      Real Field with 10 bits of precision
      sage: RealField()
      Real Field with 53 bits of precision
      sage: RealField(100000)
      Real Field with 100000 bits of precision

   Here we show the effect of rounding:

      sage: R17d = RealField(17,rnd='RNDD')
      sage: a = R17d(1)/R17d(3); a.exact_rational()
      87381/262144
      sage: R17u = RealField(17,rnd='RNDU')
      sage: a = R17u(1)/R17u(3); a.exact_rational()
      43691/131072

   Note:

     The default precision is 53, since according to the MPFR manual:
     'mpfr should be able to exactly reproduce all computations with
     double-precision machine floating-point numbers (double type in
     C), except the default exponent range is much wider and subnormal
     numbers are not implemented.'

   See also:

     * "sage.rings.real_mpfr"

     * "sage.rings.real_arb.RealBallField" (real numbers with rigorous
       error bounds)
Init docstring: Initialize self.  See help(type(self)) for accurate signature.
File:           ~/Git/sagemath/sage/src/sage/rings/real_mpfr.pyx
Type:           builtin_function_or_method

Here we create two real fields with 10 bits of precision and different rounding:

F_down = RealField(10, rnd='RNDD') # Round down
F_up = RealField(10, rnd='RNDU')   # Round up

Here we get a lower bound for $\pi$:

x_down = F_down(pi)
x_down

3.1

And an upper bound:

x_up = F_up(pi)
x_up

3.2

sme = x_down.sign_mantissa_exponent()
sme

(1, 804, -8)

xx_down = 804*2^-8
xx_down

201/64

x_up = F_up(pi)
x_up.sign_mantissa_exponent()

(1, 805, -8)

The above shows that $$804 \cdot 2^{-8} < \pi < 805 \cdot 2^{-8}.$$

We can do exact arithmetic with rationals. So, suppose we wonder if $\pi$ is bigger or smaller than $\frac{22}{7}$. We can check:

22/7 < 804*2^-8

False

805*2^-8 < 22/7

False

Okay so 22/7 was in the interval found above. (So $22/7$ is accurate to at least $2^{-8}$). We can do the same at higher precision.

To make it easier, we write afunction to convert a floating point number to the rational equal to it:

def float_to_exact_rational(x):
    s,m,e = x.sign_mantissa_exponent()
    return s*m*2^e

For example x_up equals $805 \cdot 2^{-8}$, and we can see this quickly now:

float_to_exact_rational(x_up)

805/256

From the above if we want to decide if $22/7$ is greater than or less than $\pi$, we need to compute upper and lower bounds of $\pi$ with greater precision. Let's try 20 bits of precision.

F_down = RealField(20, rnd='RNDD')
lower_bound = float_to_exact_rational(F_down(pi))
lower_bound

823549/262144

22/7 < lower_bound

False

F_up = RealField(20, rnd='RNDU')
up_bound = float_to_exact_rational(F_up(pi))
up_bound

411775/131072

up_bound < 22/7

True

We have that $\pi$ is less than up_bound (because of the rounding mode) used, and we can rigorously check that up_bound is less than $22/7$, so we have that $\pi < 22/7$.

This is an example of doing Interval Arithmetic.

Interval Arithmetic

In interval arithmetic, we represent a number $x$ by an interval $[a,b]$ that contains it. (In case we can represent a number exactly, we can represent the number as an interval consisting of a single value.)

The field RIF consists of numbers represented by intervals where the endpoints are from RR which has 53 bits of precision.

RIF

Real Interval Field with 53 bits of precision

We can construct an interval field where the endpoints have greater precision. For example:

F = RealIntervalField(100)
F

Real Interval Field with 100 bits of precision

You can construct elements in the field by converting a number into it. For example:

x = F(pi)
x

3.14159265358979323846264338328?

Numbers are printed with a question mark at the end. This question mark roughly represents the first digit which is ambiguous.

You can also construct an interval explicitly by specifing the upper and lower bounds.

x = F(2.73, 2.74)
x

2.8?

I'm not sure why the above printed 2.8?. I would think it should be 2.7?. But you can access the actual endpoints of the interval with the .lower() and .upper() methods. And these make more sense:

x.lower()

2.7299999999999999822364316059

x.upper()

2.7400000000000002131628207281

The lower bound is slightly less than $2.73$ because we round lower endpoints down, and the upper bound is slightly more than $2.74$ because we round upper endpoints up.

Here we demonstrate this using 10 bits of precision:

F = RealIntervalField(10)
F

Real Interval Field with 10 bits of precision

x = F(pi)
x

3.15?

a = float_to_exact_rational(x.lower())
a

201/64

b = float_to_exact_rational(x.upper())
b

805/256

We use x2 to store $\pi^2$:

x2 = x*x
x2

9.9?

To see the rounding moves outward observe that:

float_to_exact_rational(x2.lower()) < a^2

True

b^2 < float_to_exact_rational(x2.upper())

True

Comparisons

Recall that an equation or an inequality is only true if it is provably true. In interval arithmetic, $[a,b]$ represents any of the points in $[a, b]$. So, assuming $[a,b]$ has more than one point, $[a,b] = [a,b]$ will be false.

x2 == x2

False

If you want to check that the intervals are the same, you have to compare the lower and upper endpoints.

This can lead to some surprising situations. For example, all three of these should not be false (in normal mathematics):

9 + 18/20 < x2

False

x2 < 9 + 18/20

False

x2 == 9 + 18/20

False

Lesson: Interval fields are good (but not perfect) for testing inequalities, but not equalities! But it is really great for proving a strict inequality is true; you just need to choose a high enough precision.

Complex Analogs

A complex number is (typically represented) as $x+iy$ where $x, y \in \mathbb R$. There are analogs of the real floating point fields above, that the components in the corresponding real field.

CDF

Complex Double Field

CIF

Complex Interval Field with 53 bits of precision

ComplexField(100)

Complex Field with 100 bits of precision

F= ComplexIntervalField(100)

z = F(e^I)
z

-0.005491493961746885304119170444? - 0.042863576912993230885329718045?*I

The real and imaginary components of these complex numbers are elements of the inexact real fields discussed above. For example:

z.real().parent()

Real Interval Field with 100 bits of precision

z.imag().parent()

Real Interval Field with 100 bits of precision